#############################################################################################################################
# Illustration of MrP and MrsP with income
# November 2016
# lleemann@gmail.com
##############################################################################################################################

rm(list=ls())

# read in packages
library(foreign)
library(lme4)
library(blme)
library(arm)

# set the working directory
setwd(".../Final replication file")

# official data (referendums)
true.yesShare <- read.csv("BfS data.csv")
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# response models ---- 3 individ. predictor: GENDER, EDUCATION, AGE
model.2 <- y ~ x2.3 + x2.4 + x2.6 + x2.7 + x2.8 + x2.10 + (1|x1.gender) + (1|x1.education) + (1|x1.age) + (1|cantonnr) + (1|region)
model.L2 <- y ~ x2.3 + x2.4  + (1|x1.education) + (1|x1.gender) + (1|x1.age) + (1|cantonnr) + (1|region) + (1|x1.money)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# surevy data
voxit.data <- read.dta("Voxit.dta", convert.factors=FALSE)

vote.ids <- c(4840, 5310, 4420, 4651, 4652, 4660, 4810, 5430)


prediction.mrp <- array(NA, c(2, 26, length(vote.ids)))
prediction.mrsp <- array(NA, c(2, 26, length(vote.ids)))
N.sample.dis <- rep(NA, length(vote.ids))
pred.dis <- matrix(NA, 26, length(vote.ids))


block <- read.csv("census final BfS 11 2012 2.csv", header = TRUE)


assign("last.warning", NULL, envir = baseenv())
start.time <- Sys.time()
for (z in vote.ids){
	data.run <- voxit.data[voxit.data$AbstNrBFS==z,]
	print(z)
	official.results <- as.numeric(true.yesShare[true.yesShare[,27]==z,1:26])
	prediction.mrp[2, ,which(vote.ids==z)] <- official.results
	prediction.mrsp[2, ,which(vote.ids==z)] <- official.results

	raw.out <- data.run[,c(1,18)] 	# getting survey yes-vote and canton
	raw.out <- na.omit(raw.out)		# dropping observations which did not respond whether they voted yes or not
	N.sample.dis[which(vote.ids==z)] <- dim(raw.out)[1]
	for (p in 1:26){
		nobs.p <- length(which(raw.out[,2]==p))		# number of obs in canton p
		www <- mean(raw.out[raw.out[,2]==p,1])
		ifelse(nobs.p!=0, pred.dis[p,which(vote.ids==z)] <- www, pred.dis[p,which(vote.ids==z)] <- mean(raw.out[,1]))
		}
	data.run1 <- data.run#[,c(1,2,3,)] #[,-c(32:75)]   #372
	#data.run1 <- na.omit(data.run1)
	y <- data.run1$decx
	cantonnr <- data.run1$cantonnr
	region <- data.run1$region
	x2.3 <- data.run1$german_share
	x2.4 <- data.run1$romkath_share			
	x2.5 <- data.run1$leftvote_share
	# money
	x2.6 <- data.run1$inc_under3k_share			
	x2.7<- data.run1$inc_3kto5k_share
	x2.8 <- data.run1$inc_5kto7k_share 
	x2.9 <- data.run1$inc_7kto9k_share
	x2.10 <- data.run1$inc_over9k_share 
	# demographics
	x1.gender <- data.run1$sexe    # 1=man, 2=woman
	x1.education <- data.run1$educ
	x1.age <- data.run1$age_group
	x1.money <- data.run1$money

	summary(x1.money)

	DATT <- data.frame(x2.3,x2.4,x2.5,x2.6,x2.7,x2.8,x2.9,x2.10,cantonnr,region,y,x1.gender,x1.education,x1.age,x1.money)
	head(DATT)

		# MrP with Income shares on Level 2
		model.r <- glmer(model.2, family=binomial(probit), data=DATT, control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))
		source("source_MrP_CH_money.R")  
		prediction.mrp[1,,which(vote.ids==z)] <- prediction

		# MrsP with Income 
		model.r <- glmer(model.L2, family=binomial(probit), data=DATT, control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))
		source("source_MrsP_CH_money.R")  
		prediction.mrsp[1,,which(vote.ids==z)] <- predictionMrsP
}


end.time <- Sys.time()
cat(c("Job started at:",start.time))
cat(c("Job finished at:",end.time))
dur <- end.time-start.time
cat(paste("Job took:",dur))


diff.mrp <- prediction.mrp[1,,] - prediction.mrp[2,,]/100
diff.mrsp <- prediction.mrsp[1,,] - prediction.mrsp[2,,]/100

colMeans(abs(diff.mrp))
colMeans(abs(diff.mrsp))

pre <- abs(colMeans(abs(diff.mrp))-colMeans(abs(diff.mrsp)))/colMeans(abs(diff.mrp))
vec <- order(pre)

Title <- c("Capital Gains Tax", "New Business \n Tax Scheme", "Mileage-Based \n Heavy-Vehicle Levy", "Solar Power Tax", "Renewable \n Energy Tax", "Energy Consumption \n Tax", "Tax Energy \n Instead of Income", "Increased VAT for \n Disability Insurance" ) #"Tax Cuts for Home Owners (556)", )

## Final Plot
png("Figure05.png", width = 870, height = 700, points=16)
#par(family="CMU Serif", oma=c(0,15,0,1),mai=c(.1,.5,.1,.5),lwd=2)
par(oma=c(0,15,0,1),mai=c(.1,.5,.1,.5),lwd=2)
barplot(pre[vec], horiz=TRUE, col=rgb(0,0,255,100, maxColorValue=255), border=rgb(0,0,255,250, maxColorValue=255),xaxt="n",xlim=c(0,.8), lwd=3)
text(rep(-0.41,8),seq(0.7,9.2,length.out=8),Title,xpd=NA,pos=4, cex=1.5)
text(pre[vec]+0.05,seq(0.7,9.2,length.out=8),paste(100*round(pre[vec],2),"%"),xpd=NA, cex=1.5)
dev.off()






